Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c’est : MERVEILLEUX même si on s’arrache parfois les cheveux .
N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
# Création des différents répertoires
# Répertoire principal
mkdir -p 01_Mini_projet_Bacillus
# Direction le répertoire principal
cd ~/01_Mini_projet_Bacillus
# Création de sous-répertoires
mkdir 01_Data 02_Analysis 03_Literature 04_Archive
cd ~/01_Mini_projet_Bacillus/01_Data
mkdir -p 01_FASTQ 02_Genome 03_FASTQC 04_CLEANING 05_MAPPING
cd ~/01_Mini_projet_Bacillus/02_Analysis
mkdir -p 01_Script 02_Output
Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
## Chargement du run SRR10390685 ##
# Direction dans le répertoire ~/01_Data/01_FASTQ
cd ~/01_Mini_projet_Bacillus/01_Data/01_FASTQ
# Chargement du module sra-tools et la subcommand "fasterq-dump"
module load sra-tools
fasterq-dump
# option --split-files car séquençage paired-end
srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir 01_FASTQ
# Compression fichiers
srun gzip *.fastq
Combien de reads sont présents dans les fichiers R1 et R2 ?
## Comptage du nombre de reads en R1 et R2 ##
# Direction dans le répertoire contenant les FASTQ : ~/01_Data/01_FASTQ
cd ~/01_Mini_projet_Bacillus/01_Data/01_FASTQ
# Chargement du module seqkit
module load seqkit
# Stats sur les fichiers fastq
srun seqkit stats --threads 1 *.fastq
Les fichiers FASTQ contiennent 7,066,055 reads chacun.
Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
## Téléchargement du génome Bacillus subtilis souche ASM904v1 ##
# Dossier d'accueil
cd ~/01_Mini_projet_Bacillus/01_Data/02_Genome
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
# Decompression du fichier fasta
gunzip GCF_000009045.1_ASM904v1_genomic.fna.gz
Quelle est la taille de ce génome ?
seqkit stats --threads 1 GCF_000009045.1_ASM904v1_genomic.fna
La taille de ce génome est de 4,215,606/2 = 2,107,803 paires de bases.
Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
## Téléchargement du fichier d' annotation(.gff) Bacillus subtilis souche ASM904v1 ##
# Dossier d'accueil
cd ~/01_Mini_projet_Bacillus/01_Data/02_Genome
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
# Decompression du fichier gff
gunzip GCF_000009045.1_ASM904v1_genomic.gff.gz
Combien de gènes sont connus pour ce génome ?
grep -c "ID=gene" GCF_000009045.1_ASM904v1_genomic.gff
4536 gènes sont recensés dans le fichier d’annotation.
Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
cd ~/01_Mini_projet_Bacillus/01_Data/01_Data/01_FASTQ
# Chargement module fastqc
module load fastqc
# Run fastqc
srun --cpus-per-task 8 fastqc *.fastq.gz -o ~/01_Mini_projet_Bacillus/01_Data/03_FASTQC/ -t 8
La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?
car la qualité moyenne par base est > 30 comme le montre le rapport fastQC
Lien vers le rapport FastQC rapport FastQC
Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?
car la taille des reads n’est pas homogène sur l’ensemble du jeu de données pour R1: 130-151 et pour R2 : 35-151
Normalement, si les reads n’avaient subi aucune étape de nettoyage, il n’y aurait qu’une seul taille unique (celle du séquençage)
Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
profondeur de séquençage = nombre de bases / taille génome n
Ici : (7066055 *2)/4215606
La profondeur de séquençage est de : 3,35X.
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.
# Direction dans répertoire
cd ~/01_Mini_projet_Bacillus/01_Data
# Chargement du module fastp
module load fastp
# Run fastp
srun --cpus-per-task 8 fastp --in1 01_FASTQ/SRR10390685_1.fastq.gz --in2 01_FASTQ/SRR10390685_2.fastq.gz --out1 04_CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 04_CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html 04_CLEANING/fastp.html --thread 8 --trim_poly_g --cut_mean_quality 30 --cut_window_size 8 --trim_tail1 1 --length_required 35 --json /dev/null &> 04_CLEANING/fastp.log
Les paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
| --trim_poly_g | par défaut | car dans le rapport fastQC il y avait des reads contenant des polyG (chez Illumina, ça arrive s’il n’y a pas de signal) |
| --cut_mean_quality | 30 | tri sur la qualité moyenne des reads>=30 (standard) |
| --cut_window_size | 8 | fenêtre coulissante de 8 (3’->5’) |
| --trim_tail1 | 1 | pour supprimer la dernière base issue du dernier cycle chez Illumina, souvent de mauvaise qualité |
| --length_required | 35 | taille minimale de 35pb (choix arbitraire) |
Ces paramètres ont permis de conserver 6898098 reads pairés, soit une perte de 2,37% des reads bruts.
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].
## Alignement sur le génome de référence avec bwa ##
cd ~/01_Data/05_MAPPING
# Chargement de bwa et bwa index
module load bwa
bwa
# Indexage du génome avec bwa
bwa index
srun bwa index 02_Genome/GCF_000009045.1_ASM904v1_genomic.fna
# Run mapping
bwa mem
srun --cpus-per-task=32 bwa mem 02_Genome/GCF_000009045.1_ASM904v1_genomic.fna 04_CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz 04_CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > 05_MAPPING/SRR10390685_bwa_alignement.sam
# Chargement du module samtools
module load samtools
# Conversion du fichier sam en bam avec samtools
srun --cpus-per-task=8 samtools view --threads 8 05_MAPPING/SRR10390685_bwa_alignement.sam -b >05_MAPPING/SRR10390685_bwa_alignement.bam
# Tri du fichier bam
srun samtools sort SRR10390685_bwa_alignement.bam -o SRR10390685_bwa_alignement.sort.bam
# Indexage du fichier bam trié
srun samtools index SRR10390685_bwa_alignement.sort.bam
Combien de reads ne sont pas mappés ?
samtools view -c -f0x4 SRR10390685_bwa_alignement.sam
747443 reads ne sont pas mappés.
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
# Direction répertoire
cd ~/01_Mini_projet_Bacillus/01_Data/
# Chargement du module bedtools
module load bedtools
# Recherche du gène trmNF dans le fichier d'annotation .gff
grep trmNF GCF_000009045.1_ASM904v1_genomic.gff | awk '$3=="gene"' > trmNF_gene.gff
# Récupérez les alignements sur le gène trmNF, avec l'option -f 0.5
srun bedtools intersect -a 05_MAPPING/SRR10390685_bwa_alignement.sort.bam -b 02_Genome/trmNF_gene.gff -sorted -f 0.5 > 05_MAPPING/SRR10390685_trmNF.bam
# Tri et indexage des reads
srun samtools sort 05_MAPPING/SRR10390685_trmNF.bam -o 05_MAPPING/SRR10390685_trmNF.sort.bam
# Statistiques d'alignement avec idxstats et flagstats
srun samtools idxstats SRR10390685_trmNF.sort.bam > SRR10390685_trmNF.sort.bam.idxstats
srun samtools flagstats SRR10390685_trmNF.sort.bam > SRR10390685_trmNF.sort.bam.flagstats
2841 reads chevauchent le gène d’intérêt.
sra-tools/2.10.0
fasterq-dump version 2.10.3
seqkit v0.14.0
FastQC v0.11.9
bedtools v2.29.2
samtools 1.10, using htslib 1.10.2
bwa Version: 0.7.17-r1188
1. toolkit NS. NCBI sra toolkit. NCBI, GitHub repository. 2019.
2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Zhou Y, Chen Y, Chen S, Gu J. Fastp: An ultra-fast all-in-one fastq preprocessor. Bioinformatics. 2018;34:i884–90. doi:10.1093/bioinformatics/bty560.
4. Li H. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprint arXiv:13033997. 2013.
5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009;25:2078–9.
6. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
7. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (igv): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.
DUBii evaluation
Auteurs du TP :
Valentin Loux & Olivier Rué
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France